* load the replication data set 
* This requires the installation of both MPlus and runmplus
* MPlus needs to be in the path 
log using islamophobia-replication.txt, replace text
use islamophobia-replication.dta, replace

* run the CFA using MPlus in all five countries/regions and store results on disk 

forvalues g = 1(1)5 {
          display "group: `g'"

          runmplus pop1-pop4 cult1 cult2 econ1 econ2 auth1 auth2 auth3 islam1 islam2 belonging church if group == `g', model(POP by pop1* pop2 pop3 pop4; NAT by cult1* cult2 econ1* econ2; AUTH by auth1* auth2 auth3; RELIG by belonging* church; ISL by islam1* islam2; POP@1; NAT@1; AUTH@1; RELIG@1; ISL@1) savelog(`g'.log) savedata(save=fscores; file=fscores-`g'.dat) est(mlr)

* Extract estimates and standard errors

          matrix estimates = r(estimate)
          matrix estimates = estimates[1..56,1]
* Need to transpose this
          matrix b = estimates' 
          matrix ses = r(se)
          matrix ses = ses[1..56,1]         
* We need to square the entries in se to get back to covariances. We also remove the off-diagonal elements 
          matrix V = ses * ses'
          matrix V = diag(vecdiag(V))
* Now post it 
          ereturn post b V

          estadd scalar RMSEA = real(r(RMSEA) )
          estadd scalar TLI = real(r(TLI))
          estadd scalar CFI = real(r(CFI))
          estadd scalar N = r(Number_of_observations)   

          eststo single`g'
          }

* Now make some graphs to assess the linearity of the relationship between religion and Islamophobia/Nativism/Authoritarianism
* look only at random 10% sample to reduce overplotting          

forvalues g = 1(1)5 {
          capture frame drop fscore`g'
          frame create fscore`g'
          frame change fscore`g' 
          runmplus_load_savedata ,out (`g'.log.out) clear
          keep auth nat isl relig 
          rename isl group`g'1
          rename nat group`g'2
          rename auth group`g'3
          rename relig group`g'5
          gen byte insample = runiform()>.9    
          }


forvalues g = 1(1)1 {

frame change fscore`g'

graph twoway scatter group`g'1 group`g'5 if insample, jitter(1.3) mcolor(%30) || lowess group`g'1 group`g'5,yscale(range(-3 3)) xscale(range(-3 3)) aspect(1) xlabel(-3 0 3) ylabel(-3 0 3) ytitle(Islamophobia) xtitle(Religion) xscale(off) plotregion(style(none))  leg(off)
graph rename subgraph`g'isl , replace 

graph twoway scatter group`g'2 group`g'5 if insample, jitter(1.3) mcolor(%30) || lowess group`g'2 group`g'5,yscale(range(-3 3)) xscale(range(-3 3)) aspect(1) xlabel(-3 0 3) ylabel(-3 0 3) ytitle(Nativism) xtitle(Religion)  xscale(off) plotregion(style(none))  leg(off)
graph rename subgraph`g'nat , replace 

graph twoway scatter group`g'3 group`g'5 if insample, jitter(1.3) mcolor(%30) || lowess group`g'3 group`g'5,yscale(range(-3 3)) xscale(range(-3 3)) aspect(1) xlabel(-3 0 3) ylabel(-3 0 3) ytitle(Authoritarianism) xtitle(" ") plotregion(style(none))  leg(off)
graph rename subgraph`g'auth , replace 


}


forvalues g = 2(1)5 {

frame change fscore`g'

graph twoway scatter group`g'1 group`g'5 if insample, jitter(1.3) mcolor(%30) || lowess group`g'1 group`g'5,yscale(range(-3 3)) xscale(range(-3 3)) aspect(1) xlabel(-3 0 3) ylabel(-3 0 3) ytitle(Islamophobia) xtitle(Religion) xscale(off) yscale(off) plotregion(style(none))  leg(off)
graph rename subgraph`g'isl , replace 

graph twoway scatter group`g'2 group`g'5 if insample, jitter(1.3) mcolor(%30) || lowess group`g'2 group`g'5,yscale(range(-3 3)) xscale(range(-3 3)) aspect(1) xlabel(-3 0 3) ylabel(-3 0 3) ytitle(Nativism) xtitle(Religion) xscale(off) yscale(off) plotregion(style(none))  leg(off)  
graph rename subgraph`g'nat , replace 

graph twoway scatter group`g'3 group`g'5 if insample, jitter(1.3) mcolor(%30) || lowess group`g'3 group`g'5,yscale(range(-3 3)) xscale(range(-3 3)) aspect(1) xlabel(-3 0 3) ylabel(-3 0 3) ytitle(Authoritarianism) xtitle(" ") yscale(off) plotregion(style(none))  leg(off)
graph rename subgraph`g'auth , replace 


}

* Combine monster graph 
* This works, but if want column titles, I need to assemble the graph column-wise

graph combine subgraph1isl subgraph1nat subgraph1auth, col(1) colfirst xcommon ycommon t1title(Germany (W))
graph rename germanyw, replace 

graph combine subgraph2isl subgraph2nat subgraph2auth, col(1) colfirst xcommon ycommon t1title(Germany (E))
graph rename germanye, replace 

graph combine subgraph3isl subgraph3nat subgraph3auth, col(1) colfirst xcommon ycommon t1title(Britain)
graph rename britain, replace 

graph combine subgraph4isl subgraph4nat subgraph4auth, col(1) colfirst xcommon ycommon t1title(France)
graph rename france, replace 

graph combine subgraph5isl subgraph5nat subgraph5auth, col(1) colfirst xcommon ycommon t1title(Netherlands)
graph rename netherlands, replace 


graph combine germanyw germanye britain france netherlands ,col(5) colfirst xcommon ycommon b1title(Religion)

graph save monster, replace

graph export nonlinearities.pdf , replace 

frame change default

* Now make the tables 

esttab single1 single2 single3 single4 single5 using cov-table.tex, keep(*with*) se mlabel("Germany (W)" "Germany (E)" "Britain" "France" "Netherlands") nonum noabbrev replace b(2) stats(N RMSEA TLI CFI,fmt(a2)) varlabels(nat_with_pop "1: cov(Nativism,Populism)" auth_with_pop "2: cov(Authoritarianism,Populism)" auth_with_nat "3: cov(Authoritarianism,Nativism)" relig_with_pop "4: cov(Religiosity,Populism)" relig_with_nat "5: cov(Religiosity,Nativism)" relig_with_auth "6: cov(Religiosity,Authoritarianism)" isl_with_pop "7: cov(Islamophobia,Populism)" isl_with_nat "8: cov(Islamophobia,Nativism)" isl_with_auth "9: cov(Islamophobia,Authoritarianism)" isl_with_relig "10: cov(Islamophobia,Religiosity)") 

* Factor loadings for the appendix 
esttab single1 single2 single3 single4 single5 using loadings-table-1.tex , keep(*by_auth* *by_cult* *by_econ*) se mlabel("Germany (W)" "Germany (E)" "Britain" "France" "Netherlands") nonum noabbrev replace b(2) stats(N) subst(nat_by_ "Nativism: " auth_by_ "Authoritarianism: ")

* Second part 
esttab single1 single2 single3 single4 single5 using loadings-table-2.tex , keep(*by_pop* *by_isl*) se mlabel("Germany (W)" "Germany (E)" "Britain" "France" "NL") nonum noabbrev replace b(2) stats(N) subst(pop_by_ "Populism: " isl_by_ "Islamophobia: ")

* Third
esttab single1 single2 single3 single4 single5 using loadings-table-3.tex , keep(*by_bel* *by_ch*) se mlabel("Germany (W)" "Germany (E)" "Britain" "France" "NL") nonum noabbrev replace b(2) stats(N) subst(relig_by_ "Religiosity: ")


log close 
